#define TIME 0 //not working right now


#define KS	0
#define TSS	1


#define METRIC	TSS

#define PI	3.141592653589793238462
#define THRESHOLD	80
#define ACCEPTANCES	1000
#define REJECTIONS	ACCEPTANCES*100 //if we start rejecting 99%, consider it done.
#define GENERATIONS	33
#define PROPORTION	0.10	
#define STDEV	10
#define SFMT_MEXP 19937
#define LINE_LEN 1024// length of input lines..
#define RANDOM_SEED 1 //be careful: without this, you may run on repeat after a while!
#define INDEPENDENT 0
#define MAX(a,b)	((a)>(b) ? (a):(b))
#define MIN(a,b)	((a)<(b) ? (a):(b))
#define DELIM " \t,"
using namespace std;

#include <stdlib.h>
#include <fstream>
#include <iostream>
#include <sstream>
#include <time.h>
#include <getopt.h>
#include <math.h>
#include <string.h>
#include <assert.h>
#include <mpi.h>
#include <float.h>
#include "SFMT.h"
#include "SFMT.c"
//#include "mergesort.cpp"
#include "rand.cpp"
#include "Rmn.cpp"
#include "mpi.cpp"
#include "mergesort.cpp"
#include "customizable.cpp"


struct Rmn *generate_emperical(uint n, uint m) {
struct Rn *true_mean=init_Rn(n,1),*mean=NULL; 
struct Rmn *forced_covariance=init_Rmn(n,n,0),*covariance=NULL;

for (uint i=0;i<n;i++)
	for (uint j=0;j<=i;j++) {
		forced_covariance->f[i]->f[j]=0.75+0.1*u01();
		forced_covariance->f[j]->f[i]=forced_covariance->f[i]->f[j];;
	}

struct svd *true_ed=svd_svd(forced_covariance);//=new svd;

struct Rmn *d_emperical=theorize_Rmn(true_mean,true_ed,m), *true_covariance=theorize_covariance_Rmn(true_ed);


	cout<<"#True mean vector is:"<<endl;
	write_Rn_comment(true_mean,cout);
	cout<<"#True eigen values are:"<<endl;
	write_Rn_comment(true_ed->values,cout);
	cout<<"#True eigenvectors are:"<<endl;
	write_Rmn_comment(true_ed->vectors,cout);
	cout<<"#True covariance matrix is:"<<endl;
	write_Rmn_comment(true_covariance,cout);
	cout<<"#True r^2 matrix is:"<<endl;
	write_Rmn_comment(r2_Rmn(true_covariance),cout);
	cout<<"#Inverse of true covariance matrix is:"<<endl;
	write_Rmn_comment(theorize_inverse_covariance_Rmn(true_ed),cout);
	cout<<"#And multiplied those are (sanity check...)"<<endl;
	write_Rmn_comment(multiply_Rmn(true_covariance,theorize_inverse_covariance_Rmn(true_ed)),cout);

	cout<<"#Has cholesky decomp:"<<endl;
	write_Rmn_comment(cholesky_Rmn(true_covariance),cout);

	struct svd *tmp_eigen=svd_svd(true_covariance);
	cout<<"#Has SVD matrix:"<<endl;
	write_Rmn_comment(tmp_eigen->vectors,cout);
	cout<<"#And SVD vector:"<<endl;
	write_Rn_comment(tmp_eigen->values,cout);
	cout<<"#Putting those two back together, we have:"<<endl;
	write_Rmn_comment(theorize_covariance_Rmn(tmp_eigen),cout);
	cout<<"#with theoretical determinant: "<<theorize_determinant(true_ed->values)<<endl;
	cout<<"#and true determinant:"<<calculate_determinant(true_covariance)<<endl;
	//GENERATE THE EMPERICAL DATASET HERE...
	cout<<"#prior prob of mu at peak of integral is:"<<mvn_pdf(true_mean,true_ed,true_mean)<<endl;
	cout<<"#prior prob of s at peak of integral is:"<<inverse_wishart_pdf(true_ed,true_ed)<<endl;

	 mean=mean_Rn(d_emperical); //get the mean of the emperical dataset
        covariance=cov_Rmn(d_emperical);

	cout<<"#Emperical dataset has mean vector:"<<endl;
	write_Rn_comment(mean,cout);
	cout<<"#and covariance matrix:"<<endl;
	write_Rmn_comment(covariance,cout);
	cout<<"#and r2 matrix:"<<endl;
	write_Rmn_comment(r2_Rmn(covariance),cout);

return d_emperical;

}


int main(int argc, char *argv[]) { 

/*
float random_data[10];
for (uint i=0;i<10;i++)
		random_data[i]=u01();
cerr<<"here's the random data:"<<endl;
for (uint i=0;i<10;i++)
		cerr<<random_data[i]<<endl;
uint *index=mergeSort(random_data,10);
cerr<<"here's the sorted data:"<<endl;
for (uint i=0;i<10;i++)
		cerr<<random_data[i]<<"\t"<<index[i]<<endl;
*/
	//cerr<<"ACCEPTANCES/PROPORTION="<<ACCEPTANCES/PROPORTION<<endl;
	uint m=20,n=1,random;
int optc;

struct Rmn *d_emperical=NULL;

static struct option const longopts[] = {
                {"device", required_argument, NULL, 'd'},
		{"generate",required_argument,NULL,'g'},
                {"help", no_argument, NULL, 'h'},
                {"version",no_argument,NULL,'v'},
                {NULL,no_argument,NULL,0} };

while ((optc=getopt_long (argc, argv, "d:g:hv", longopts, NULL)) != -1)
        switch (optc) {
                case 'd': //select device type... for future OpenCL implementation
                        if (strcasecmp("cpu",optarg)==0) {
                                fprintf(stderr,"Setting cl_device_type to CL_DEVICE_TYPE_CPU...\n");
                             //   DEVICE_TYPE=CL_DEVICE_TYPE_CPU;
                        } else if (strcasecmp("gpu",optarg)==0) {
                                fprintf(stderr,"Setting cl_device_type to CL_DEVICE_TYPE_GPU...\n");
                             //   DEVICE_TYPE=CL_DEVICE_TYPE_GPU;
                        } else if (strcasecmp("all",optarg)==0) {
                                fprintf(stderr,"Setting cl_device_type to CL_DEVICE_TYPE_ALL...\n");
                             //   DEVICE_TYPE=CL_DEVICE_TYPE_ALL;
                        } else {
                                fprintf(stderr,"Unrecognized cl_device type '%s'. Leaving as CL_DEVICE_TYPE_GPU...\n",optarg);
                        }
                        break;
		case 'g': //generate test data and exit
						n=atoi(strtok(optarg,","));
			m=atoi(strtok(NULL,","));

			assert(n>0 && m>0);
			if (mpi_rank==0) {
				d_emperical=generate_emperical(m,n);
				write_Rmn(d_emperical,cout);
			}
			exit(EXIT_SUCCESS);
			break;
		case 'h': // print help menu and exit;
			cerr<<"Help!"<<endl;
			exit(EXIT_SUCCESS);
			break; 
		case 'v': 
			cerr<<"Version 2013-10"<<endl; exit(EXIT_SUCCESS);
			break;
	}



MPI::Init();
bool quit=0;
mpi_rank=MPI::COMM_WORLD.Get_rank(), mpi_ranks=MPI::COMM_WORLD.Get_size(), accepts_per_rank=ACCEPTANCES/mpi_ranks;
assert(ACCEPTANCES%mpi_ranks==0);
float rejection_threshold=0,d, max_weight=1.f/ACCEPTANCES,min_weight=1.f/ACCEPTANCES,sum_weight=1, sorted_distances[(uint)(ACCEPTANCES/PROPORTION)];//,last_weight_var=0;

uint *sorted_distances_index=new uint[ACCEPTANCES];
for (uint i=0;i<ACCEPTANCES;i++)
		sorted_distances_index[i]=i;

//true_ed->values=init_eigen_Rn(n);
//true_ed->vectors=orthogonalize_Rmn(init_Rmn(n,n,1));

//ERASE...
//for (uint i=0;i<n;i++)
//	for (uint j=0;j<n;j++)
	//		if (i==j)
//			true_ed->values->f[i]=1;//,true_ed->vectors->f[i]->f[j]=1;
	//	else
	//		true_ed->vectors->f[i]->f[j]=1+u01();


struct Rn *tmp_Rn=NULL, *mean=NULL, *sample_mean=NULL, *last_ecdf=init_Rn(ACCEPTANCES,0), *ecdf=init_Rn(ACCEPTANCES,0);
struct Rmn *tmp_Rmn=NULL, *d_theoretical=NULL, *covariance=NULL;

struct mpi_progress mpi_progress=init_mpi_progress(n,ACCEPTANCES/PROPORTION);

       struct svd *sample_ed=NULL;
if (d_emperical==NULL)  {
	char line[LINE_LEN], *ptr_a=NULL;
	float data[LINE_LEN][LINE_LEN];

	for  (m=0;fgets(line,LINE_LEN,stdin);) {
		if (line[0]!='#') {
			data[m][0]=atof(strtok(line,DELIM)); // first 
			for (n=1;(ptr_a=strtok(NULL,DELIM))!=NULL;n++)
				data[m][n]=atof(ptr_a);
			m++;		
		}
	}

	d_emperical=init_Rmn(m,n,0);
	for (uint i=0;i<m;i++)
		for (uint j=0;j<n;j++)
			d_emperical->f[i]->f[j]=data[i][j];
}

	cerr<<"dimension of data in m="<<m<<",n="<<n<<endl;
	cerr<<"data is:"<<endl;
	write_Rmn(d_emperical,cerr);

//	cout<<"Normalized dataset is:"<<endl;
//	write_Rmn_comment(normalize_Rmn(d_emperical),cout);



if (METRIC==TSS) {
	for (uint i=0;i<m*m;i++) {
		rejection_threshold+=distance_Rmn(tmp_Rmn=permute_Rmn(cp_Rmn(d_emperical)),d_emperical);
		destroy_Rmn(tmp_Rmn);
	}
	rejection_threshold/=(float)(m*m);
	cerr<<"Initial rejection threshold is:"<<rejection_threshold<<endl;
}
else if (METRIC==KS)
	rejection_threshold=1;

 mean=mean_Rn(d_emperical); //get the mean of the emperical dataset
        covariance=cov_Rmn(d_emperical);
	
	struct svd *ed=svd_svd(covariance);

if (mpi_rank==0) {
	/*
        ostringstream s_output;
        s_output<<"analytical.txt";
        ofstream f_output;
        f_output.open(s_output.str().data(),ios::out);
        assert(f_output.is_open());
//        write_Rn(true_mean,f_output);
//	write_Rmn(true_covariance,f_output);
        f_output.close();


	//this is ghetto... go back and do it correctly
	ostringstream s_output2;
	s_output2<<"emperical.txt";
	ofstream f_output2;
	f_output2.open(s_output2.str().data(),ios::out);
	assert(f_output2.is_open());
	write_Rmn(d_emperical,f_output2);
	f_output2.close();
*/
}



struct Rmn **last_eigenvectors=new Rmn*[mpi_ranks*accepts_per_rank], *last_mean=theorize_Rmn(mean,ed,ACCEPTANCES/PROPORTION),/*init_Rmn(mpi_ranks*accepts_per_rank,n,1)*/ *last_eigenvalues=init_Rmn(mpi_ranks*accepts_per_rank,n,0);


for (uint i=0;i<mpi_ranks*accepts_per_rank;i++) {
	destroy_Rn(last_eigenvalues->f[i]);
	last_eigenvalues->f[i]=init_eigen_Rn(n);
	last_eigenvectors[i]=init_Rmn(n,n,1);
}

//uint reset_count;

for (uint T=0;T<GENERATIONS;T++) {

	//start: ;

mpi_progress.start_time=MPI::Wtime();

//	reset_count=0;

	mpi_progress.acceptances[mpi_rank]=0;
	mpi_progress.rejections[mpi_rank]=0;

	tmp_Rmn=cov_Rmn(last_mean);

	//what happens if we zero out the off diagonals?
	if (INDEPENDENT)
		for (uint i=0;i<n;i++)
			for (uint j=0;j<n;j++)
				if (i!=j)
					tmp_Rmn->f[i]->f[j]=0;
	
	scale_Rmn(tmp_Rmn,2); // paper recommended doubilng the variance...

	destroy_svd(sample_ed);
	sample_ed=svd_svd(tmp_Rmn);

	destroy_Rmn(tmp_Rmn);



uint count=0;

	while (1) {

		while ((random=sorted_distances_index[(uint)(u01()*mpi_ranks*accepts_per_rank)])) { // pick an element at random, uniformly from the sorted_distances array 
			if (u01()<=mpi_progress.last_weights[random]/max_weight) // compare the weight to another U=u01() to see if we should pick it. divide by max_weight to make this quicker... 
				break; 
		}

		destroy_Rn(sample_mean);	
		sample_mean=perturb_Rn(last_mean->f[random],sample_ed,1);

		destroy_Rmn(d_theoretical);
		d_theoretical=theorize_Rmn(sample_mean,sample_ed,m);
	
		d=distance_Rmn(d_emperical,d_theoretical);
	
		if ( d<rejection_threshold || (PROPORTION<1  && T>0 ))  {
			
			//cout<<"this dataset has prior prob: "<<mvn_pdf(sample_mean,ed,mean)<<endl;
	
count=0;	
			//give them all weights, whether they are actually used will be seen.
			mpi_progress.current_weights[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]=0;
			for (uint i=0;i<mpi_ranks*accepts_per_rank/PROPORTION;i++) 
				if (mpi_progress.accepted[i]) //this limits it to those accepted weights from last generation 
					{
					mpi_progress.current_weights[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]+=mpi_progress.last_weights[i]*mvn_pdf(sample_mean,sample_ed,mpi_progress.mean->f[i]);
					count++;
					}
			mpi_progress.current_weights[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]=mvn_pdf(sample_mean,ed,mean)/mpi_progress.current_weights[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]];
			
//cerr<<"count is:"<<count<<", accepts per rank/PROPORTION="<<accepts_per_rank/PROPORTION<<endl;	

			mpi_progress.rejection_threshold[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]=d;
			
			destroy_Rn(mpi_progress.mean->f[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]);
//			free(mpi_progress.mean->f[mpi_rank*accepts_per_rank+mpi_progress.acceptances[mpi_rank]]);
			mpi_progress.mean->f[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]=cp_Rn(sample_mean);


			destroy_Rn(mpi_progress.eigenvalues->f[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]);
//			free(mpi_progress.eigenvalues->f[mpi_rank*accepts_per_rank+mpi_progress.acceptances[mpi_rank]]);
			mpi_progress.eigenvalues->f[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]=cp_Rn(sample_ed->values);

			destroy_Rmn(mpi_progress.eigenvectors[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]);
//			free(mpi_progress.eigenvectors[mpi_rank*accepts_per_rank+mpi_progress.acceptances[mpi_rank]]);
			mpi_progress.eigenvectors[(uint)(mpi_rank*accepts_per_rank/PROPORTION)+mpi_progress.acceptances[mpi_rank]]=cp_Rmn(sample_ed->vectors);
			
			mpi_progress.acceptances[mpi_rank]++;
			if (mpi_progress.acceptances[mpi_rank]>=accepts_per_rank/PROPORTION)  {
			//without this, we have a memory leak. with it, we have a seg fault.
			//	destroy_Rn(mean); 
		//		destroy_Rmn(covariance);
		//		destroy_Rmn(d_theoretical);
				break;	
			}	
		} else  {
			mpi_progress.rejections[mpi_rank]++;
			if (REJECTIONS>0 && mpi_progress.rejections[mpi_rank]>=REJECTIONS) {
				quit=1; break;
	

                /*		rejection_threshold+=0.1f/GENERATIONS;
				cerr<<"reset rejection_threshold to:"<<rejection_threshold<<endl;
				mpi_progress.rejections[mpi_rank]=0;
				reset_count++;
				if (reset_count==10) {
					quit=1;
					break;
				}
		//		goto start; //break;
				*/
			}
			if (mpi_rank==0) 
				if (mpi_progress.rejections[mpi_rank]%10000==0) {
					cerr<<"Iteration "<<T+1<<"/"<<GENERATIONS<<" "<<100.f*(mpi_progress.acceptances[mpi_rank]/(float)accepts_per_rank)<<"% complete, ~"<<(1-mpi_progress.acceptances[mpi_rank]/(float)accepts_per_rank)/(mpi_progress.acceptances[mpi_rank]/(float)accepts_per_rank)*(MPI::Wtime()-mpi_progress.start_time)<<"s remaining\r";
//, expected to be done in "<< (MPI::Wtime()-mpi_progress.elapsed_time[mpi_rank])/((mpi_progress.acceptances[mpi_rank]/(float)accepts_per_rank))-MPI::Wtime() <<" seconds\r";
						//<<100.f*(mpi_progress.rejections[mpi_rank]/(float)(mpi_progress.rejections[mpi_rank]+ mpi_progress.acceptances[mpi_rank]))<<"% are rejections\r";
				}
		}

	}


	/*for (uint i=0;i<(uint)(acceptances/PROPORTION);i++) {
		mpi_progress.current_weights[i]=0;
		if (mpi_progress.accepted[i])
			for (j=0;j<(uint)(acceptances/PROPORTION);j++)
			mpi_progress.current_weights[i]+=mpi_progress.last_weights[j]*mvn_pdf()
	}*/

	if (PROPORTION < 1) { 

		for (uint i=0;i<(uint)(ACCEPTANCES/PROPORTION);i++) {
				sorted_distances[i]=mpi_progress.rejection_threshold[i];
				mpi_progress.accepted[i]=0;
		}
	//	if (sorted_distances_index!=NULL)
		free(sorted_distances_index);
		sorted_distances_index=mergeSort(sorted_distances,ACCEPTANCES/PROPORTION);
		for (uint i=0;i<ACCEPTANCES;i++) 
			mpi_progress.accepted[sorted_distances_index[i]]=1; 
	}



	mpi_progress.elapsed_time=MPI::Wtime();
	min_weight=1, max_weight=0,sum_weight=0;
	for (uint i=0;i<accepts_per_rank/PROPORTION;i++)
		if (mpi_progress.accepted[mpi_rank*(uint)(accepts_per_rank/PROPORTION)+i]) {
			sum_weight+=mpi_progress.current_weights[mpi_rank*(uint)(accepts_per_rank/PROPORTION)+i];

		}
	MPI::COMM_WORLD.Allreduce(&sum_weight,&sum_weight,1,MPI::FLOAT,MPI_SUM);


	for (uint i=0;i<accepts_per_rank/PROPORTION;i++)  
		if (mpi_progress.accepted[mpi_rank*(uint)(accepts_per_rank/PROPORTION)+i]) {
			mpi_progress.current_weights[mpi_rank*(uint)(accepts_per_rank/PROPORTION)+i]/=sum_weight; //this scales everything back to 1
			min_weight=MIN(min_weight,mpi_progress.current_weights[mpi_rank*(uint)(accepts_per_rank/PROPORTION)+i]);
			max_weight=MAX(max_weight,mpi_progress.current_weights[mpi_rank*(uint)(accepts_per_rank/PROPORTION)+i]);
		}
	MPI::COMM_WORLD.Allreduce(&max_weight,&max_weight,1,MPI::FLOAT,MPI_MAX); //broadcast max_weight...
	MPI::COMM_WORLD.Allreduce(&min_weight,&min_weight,1,MPI::FLOAT,MPI_MIN);
	cerr<<"min weight is: "<<min_weight<<", max weight is: "<<max_weight<<" and sum weight is"<<sum_weight<<endl;
	assert(0<=min_weight && max_weight<=1);


/*
	struct Rn *tmp_current=init_Rn(ACCEPTANCES,0);
	struct Rn *tmp_last=init_Rn(ACCEPTANCES,0);	
	delete tmp_current->f;
	delete tmp_last->f;
	tmp_current->f=new float[ACCEPTANCES];
	tmp_last->f=new float[ACCEPTANCES];
	memcpy(tmp_current->f,mpi_progress.current_weights,sizeof(float)*ACCEPTANCES);
	memcpy(tmp_last->f,mpi_progress.last_weights,sizeof(float)*ACCEPTANCES);
*/
	
	//last_weight_var=variance_f(mpi_progress.last_weights,ACCEPTANCES);


        //if (T>0)                                        // i doubled this since the wiki article seems to hint at this in the critical value table... and change sqrt(ACCEPTANCES)
                //cout<<"K-S p-value on variable 1 is:'"<<2.f*kolmogorov_p(sqrt(2.f*ACCEPTANCES/(1.f*ACCEPTANCES*ACCEPTANCES))*kolmogorov_smirnov_compare(last_ecdf,ecdf))<<"'\n"<<endl;
                //cout<<"K.S. critical value on weights is: '"<<kolmogorov_smirnov_compare(array_to_Rn(mpi_progress.last_weights,ACCEPTANCES),array_to_Rn(mpi_progress.current_weights,ACCEPTANCES))<<"'\n"<<endl;



	for (uint j=0;j<accepts_per_rank/PROPORTION;j++)   // copy the current weights to the last weights..
		if (mpi_progress.accepted[j])
			mpi_progress.last_weights[mpi_rank*accepts_per_rank+j]=mpi_progress.current_weights[mpi_rank*accepts_per_rank+j];
		else
			mpi_progress.last_weights[mpi_rank*accepts_per_rank+j]=0;

	MPI::COMM_WORLD.Allreduce(&mpi_progress.elapsed_time,&mpi_progress.elapsed_time,1,MPI::DOUBLE,MPI_MAX);

//FIXME
	 for (uint i=0;i<mpi_ranks;i++) {
		MPI::COMM_WORLD.Bcast(&mpi_progress.rejections[i],1,MPI::INT,i);
		MPI::COMM_WORLD.Bcast(&mpi_progress.acceptances[i],1,MPI::INT,i);

		for (uint j=0;j<accepts_per_rank;j++)  {
			MPI::COMM_WORLD.Bcast(&mpi_progress.rejection_threshold[i*accepts_per_rank+j],1,MPI::FLOAT,i);
			MPI::COMM_WORLD.Bcast(&mpi_progress.last_weights[i*accepts_per_rank+j],1,MPI::FLOAT,i);
			


			 MPI::COMM_WORLD.Bcast(&mpi_progress.mean->f[i*accepts_per_rank+j]->f[0],sample_mean->n,MPI::FLOAT,i);
			 MPI::COMM_WORLD.Bcast(&mpi_progress.mean->f[i*accepts_per_rank+j]->n,1,MPI::FLOAT,i);

			 MPI::COMM_WORLD.Bcast(&mpi_progress.eigenvalues->f[i*accepts_per_rank+j]->f[0],sample_ed->values->n,MPI::FLOAT,i);
			 MPI::COMM_WORLD.Bcast(&mpi_progress.eigenvalues->f[i*accepts_per_rank+j]->n,1,MPI::FLOAT,i);


			for (uint k=0;k<n;k++)
				MPI::COMM_WORLD.Bcast(&mpi_progress.eigenvectors[i*accepts_per_rank+j]->f[k]->f[0],sample_mean->n,MPI::FLOAT,i);
			MPI::COMM_WORLD.Bcast(&mpi_progress.eigenvectors[i*accepts_per_rank+j]->n,1,MPI::FLOAT,i);

		}
	}

//	if (check_mpi_progress(mpi_progress,TIME,THRESHOLD,REJECTIONS)<0) // if it trips a wire, we quit at the end of this iteration
//		quit=1;
	
//
	if (METRIC==TSS)
		rejection_threshold=new_threshold(mpi_progress);
	else if (METRIC==KS)
		rejection_threshold-=1.f/GENERATIONS;//=(1-(T+1)/(float)GENERATIONS);//*=0.95;//(1-(T+1)/(float)GENERATIONS);//new_threshold(mpi_progress);
//FIXME

	for (uint i=0;i<ACCEPTANCES;i++) {
		ecdf->f[i]=mpi_progress.mean->f[i]->f[0];
	}

//	if (T>0) 					
//		cout<<"K-S p-value on variable 1 is:'"<<kolmogorov_p(kolmogorov_smirnov_compare(last_ecdf,ecdf),ACCEPTANCES)<<"' with critical value:'"<<kolmogorov_smirnov_compare(last_ecdf,ecdf)<<"'\n"<<endl;
	destroy_Rn(last_ecdf);
	last_ecdf=cp_Rn(ecdf);


	for (uint i=0;i<mpi_ranks*accepts_per_rank;i++)  { //this is stupid... this needs to be transposed
		destroy_Rn(last_mean->f[i]);
		last_mean->f[i]=cp_Rn(mpi_progress.mean->f[i]);
		destroy_Rmn(last_eigenvectors[i]);
		last_eigenvectors[i]=cp_Rmn(mpi_progress.eigenvectors[i]);
		destroy_Rn(last_eigenvalues->f[i]);
		last_eigenvalues->f[i]=cp_Rn(mpi_progress.eigenvalues->f[i]);
	}

	if (mpi_rank==0) {

		cout<<"After "<<mpi_progress.elapsed_time-mpi_progress.start_time<<" seconds, T="<<T<<", new threshold is: "<<rejection_threshold<<endl;
		cout<<"Generated dataset has the mean vector:"<<endl;


        ostringstream s_output;
	s_output<<"simulated.txt"<<T;
        ofstream f_output;
        f_output.open(s_output.str().data(),ios::out);
        assert(f_output.is_open());
        write_Rmn(last_mean,f_output);
        f_output.close();


		write_Rn(tmp_Rn=mean_Rn(last_mean),cout);
		destroy_Rn(tmp_Rn);

		cout<<"And covariance matrix:"<<endl;
		write_Rmn(tmp_Rmn=cov_Rmn(last_mean),cout);
		destroy_Rmn(tmp_Rmn);
			
		s_output<<"weights";
		f_output.open(s_output.str().data(),ios::out);
		assert(f_output.is_open());
		for (uint i=0;i<ACCEPTANCES;i++) 
			f_output<<mpi_progress.last_weights[i]<<endl;
		f_output.close();

//        cerr<<"K-S p-value against U(0,ACCEPTANCES) is: "<<kolmogorov_p(sqrt(ACCEPTANCES)*kolmogorov_smirnov(mpi_progress.current_weights,ACCEPTANCES))<<endl;
//http://www.real-statistics.com/statistics-tables/kolmogorov-smirnov-table/ says we need 1.63/sqrt(ACCEPTANCES) for ACCEPTANCES>50 and alpha=0.01
	// http://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test says that sqrt(n)D_n>K_alpha is our test statistic. So I guess K_0.01=1.63

				//cerr<<"X^2 statistic is: "<<(ACCEPTANCES-1)*variance_f(mpi_progress.current_weights,ACCEPTANCES)/last_weight_var<<" with "<<ACCEPTANCES-1<<" degrees of freedom"<<endl;

			        //cerr<<"Max distance between current & last weights is: "<<variance_f(mpi_progress.current_weights,ACCEPTANCES)<<" and "<<last_weight_var<<endl;


//		cout<<"And mean eigen values:"<<endl;
//		write_Rn(tmp_Rn=mean_Rn(last_eigenvalues),cout);
//		destroy_Rn(tmp_Rn);
//		cout<<"and mean eigen vectors:"<<endl;
//		write_Rmn(tmp2_Rmn=orthogonalize_Rmn(tmp_Rmn=Rmn_averages(last_eigenvectors,mpi_ranks*accepts_per_rank)),cout);
//		destroy_Rmn(tmp_Rmn); destroy_Rmn(tmp2_Rmn);

		cout<<"Acceptance rate was: "<<mpi_progress.acceptances[mpi_rank]/((float)(mpi_progress.rejections[mpi_rank]+ACCEPTANCES))<<endl;
        mpi_progress.rejections[mpi_rank]=0;


		cout<<"**************************************"<<endl;
	}
	if (quit)
		break;
	}
MPI::Finalize();
exit(EXIT_SUCCESS); 

}
